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D I We study the behaviour of the Binder cumulant on finite square 

^ ■ lattices at the Kosterlitz-Thouless phase transition. We determine the 

fixed point value of the Binder cumulant and the coefficient of the lead- 
ing logarithmic correction. These calculations are supplemented with 
Monte Carlo simulations of the classical XY (plane rotator) model, 
■ the Villain model and the dual of the absolute value solid-on-solid 

^ ! model. Using the single cluster algorithm, we simulate lattices up 

to L = 4096. For the lattice sizes reached, subleading corrections 
^ ■ are needed to fit the data for the Binder cumulant. We demonstrate 

O . that the combined analysis of the Binder cumulant and the second 

moment correlation length over the lattice size allows for an accurate 
determination of the Kosterlitz-Thouless transition temperature on 
, relatively small lattices. We test the new method at the example of 

' the 2-component (p^ model on the lattice. 
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1 Introduction 



Two dimensional systems with short range interactions and 0(2) symmetry 
undergo a Kosterhtz-Thouless (KT) phase transition [Ij. This phase tran- 
sition is of particular interest because of its peculiar nature and the large 
number of effectively two dimensional systems that are supposed to undergo 
such a transition. 

Following the theorem of Mermin and Wagner [2] the magnetisation of 
two dimensional systems with short range interactions and a continuous sym- 
metry vanishes at any finite temperature. Kosterlitz and Thouless |lj have 
argued that nevertheless two dimensional systems with 0(2) symmetry un- 
dergo a phase transition at a finite temperature. In the low temperature 
phase, the decay of the correlation function follows a power law. At a suffi- 
ciently high temperature, pairs of vortices unbind and disorder the system, 
resulting in a finite correlation length ^; i.e. the correlation function decays 
exponentially. Starting from the seminal work of Kosterlitz and Thouless a 
rather solid theoretical understanding of the transition has been established. 
See e.g. refs. [Hll]. Furthermore, there are exactly solved models [5], [6l [H E] 
that display the behaviour predicted by KT-theory. 

This excellent theoretical understanding is contrasted by the fact that 
Monte Carlo studies of the KT-transition are notoriously difficult. These 
difficulties are related with logarithmic corrections that are present in the 
neighbourhood of the transition. 

In Monte Carlo studies of critical phenomena, renormalization group in- 
variant quantities, which are also called phenomenological couplings, are very 
useful tools. In finite size scaling (FSS), they allow to locate the transition 
point, to determine the nature of the transition and, in the case of a sec- 
ond order transition, to determine the critical exponent v of the correlation 
length. The prototype of such a quantity is the so called Binder cumulant 

c/=«!5^ (1) 

where m is the magnetisation of the system 0. While the Binder cumulant is 
a standard tool in the study of second order transitions, only a few authors 
have advocated its use in the case of a KT-transition [9l [10]. This is due to 
the fact that little is known about the behaviour of the Binder cumulant in 

*The standard convention is J7 = 1 — | '^'^w-fi)'^ ■ 
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the neighbourhood of the KT-transition. The main purpose of the present 
paper is to fill this gap. 

The outline of the paper is the following: First we define the models and 
the observables that we study. Next we summarize the results from KT- 
theory which are relevant to our problem. Then we derive the fixed point 
value of the Binder cumulant and the coefficient of the leading logarithmic 
correction. We discuss our Monte Carlo simulations of the XY model, the 
Villain model and the dual of the absolute value solid-on-solid (ASOS) model 
at the KT-transition. The Monte Carlo results for the Binder cumulant, the 
second moment correlation length and the helicity modulus are confronted 
with the theoretical predictions. Based on this discussion, we propose to 
determine the transition temperature of a model which is supposed to un- 
dergo a KT-transition by simultaneously matching the values of the second 
moment correlation length and the Binder cumulant with those obtained for 
the three models studied in this paper. This method is tested at the example 
of the two-component 0^ model on the square lattice. Finally we give our 
conclusions. 

2 The models 

We study the XY model on the square lattice and two generalizations of it: 
The Villain model and the dual of the absolute value solid-on-solid model. 
This allows us to check the universality of the behaviour of the Binder cumu- 
lant and other phenomenological couplings. The Boltzmann factor of these 
models can be written as a product of weights for pairs of nearest neighbour 
sites on the lattice. In all three cases we have 

Bi{s,}) = l[w{9,,,) , (2) 

where 9x,^ is the angle between nearest neighbour spins s^: and s^+p,- The 
spin Sx is a unit vector with two real components, x = {xi,X2) labels the 
sites on the square lattice, where xi G {1, 2, Li} and X2 G {1, 2, L2} Q, 
fi gives the direction on the lattice and /i is a unit-vector in the /x-direction. 
We consider periodic boundary conditions in both directions. The weight 
function is periodic: w{6) = w{6 + 2n7r) for any integer n. The partition 

^In our simulations we use Li ^ L2 — L throughout 
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function is given by the integral 

Z = [ n[d4'^d.f'^(^"'^' - 1)] B{{s.}) ■ (3) 

Note that two dimensional XY-models on square lattices with Boltzmann 
factors given by eq. ([2]) can be exactly mapped onto so called solid-on-solid 
(SOS) models. For a detailed discussion see ref. |Tl]. The variables of 
an SOS model are integers living on the sites a; of a square lattice. The 
Boltzmann factor of these models can be written as a product over nearest 
neighbour pairs 

B{{hx}) = \{w{\K-h,+ -,\) , (4) 

X,fl 

The relation between the weight function of an SOS model and its dual XY 
model is given by 

oo 

w{9) = w{n) cos{n6) (5) 

n=— oo 

and ^ 

w(n) = — [ d9 w{9) cosine) . (6) 

The weight of the standard XY model (or plane rotator model) is given 

by 

wxY{dx,p,) = exp(/3 cos(6'^,^)) = exp(/5 s- s^+f,) , (7) 

where (3 is the inverse temperature. In the case of the other two models it 
is simpler to give the weights for the SOS representation. The Villain model 
is dual to the discrete Gaussian solid-on-solid (DGSOS) model. Its weight is 
given by 

WDGSOsin) = exp(-^n2) . (8) 

It follows 

oo ^ oo . ^ V 

wviUainiO) = ^ exp{-Pn^ + in9) = c ^ exp f --(6* - 27rm)M , (9) 

11 = — rv~i m = — rv^ ^ ^ 



where P = 1/(2/3) and c = y/27ip. For a proof of this equation see e.g. 
appendix A.l of ref. [12]. 
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The weight of the absolute value solid-on-solid (ASOS) model is given by 



WAsos{n) = exp{-(3\n\) . (10) 

The KT-temperatures of these models have been accurately determined 
[13, HI] by an RG mapping of these models with the exactly solved body 
centred cubic solid-on-solid (BCSOS) model P, El El E] : 



= 1.1199(1) 6^^ = 0.93(1), 

^villa^n ^ 0.75154(23) 6^^""*" = 0.32(1) , 
fjf^os ^ 0.80608(2) bi^'^^ = 2.78(3) . (11) 

In addition to the (inverse) temperature of the KT-transition we give the 
scale factor bm that is needed to match the given model with the BCSOS 
model. E.g. results for the ASOS model at the KT-transition obtained at 
the lattice size L^sos match with the BCSOS model at the transition at the 
lattice size Lbcsos = LAsos/^m^^- 



3 The observables 

In this section we shall summarize the definitions of the observables that we 
have measured in our simulations. The total magnetisation is defined by 



m 

X 



The magnetic susceptibility is then given by 

X=^{m') , (13) 

where (. . .) denotes the expectation value with respect to the Boltzmann fac- 
tors defined in the previous section. For completeness we repeat the definition 
of the Binder cumulant 

(_(mY) 

The second moment correlation length on a lattice of the size is defined 

by 

1 / y \ 1/2 
= TT— TTTTT ( I - 1 ) . (IS) 



U-^- (14) 



2sin(7r/L) \F 
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where x is the magnetic susceptibihty as defined above and 



F = ^(44) cos(27r(yi - (16) 
x,y 

is the Fourier transform of the correlation function at the smallest non- 
vanishing momentum. In our simulations we have measured F for both 
directions of the lattice to reduce the statistical error. 

The helicity modulus T gives the reaction of the system under a torsion 
[T5] . To define the helicity modulus we consider a system, where rotated 
boundary conditions in one direction are introduced: For pairs x, x + /t of 
nearest neighbour sites on the lattice with xi = Li and /i = 1 the weight 
w{6x,^i) is replaced by w{6x,^ + a). The hehcity modulus is then defined by 
the second derivative of the free energy with respect to a at a = 



L2 da^ 



(17) 

a=0 



Note that we have skipped a factor one over temperature in our definition of 
the helicity modulus to obtain a dimensionless quantity. In the case of the 
standard XY model it is easy to write the helicity modulus as an observable 
of the system at a = [16]. For Li = L2 = L we get 

T=|(j.j.,o-f^((^'e-^«)^) . (18) 

Note that this equation is not valid for the other two models that we have 
studied. 



4 KT-theory 

At low temperatures, fiuctuations are suppressed and we might expand the 
weight as 

- ln(u'(^)) = const - + . . . (19) 

Note that for the models discussed above, w{9) is an even function that 
assumes its maximum at 6' = 0. Using this approximation we arrive at the 
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exactly solvable Gaussian model (or free field theory in the language of high 
energy physics): 



J llW6 f 5^ ^ J exp{-Hsw) (20) 



with 

Hsw = — ^^{^x — (^x+fif ■ (21) 

Note that the 5-function in eq. (I20l) is needed to render the integral finite. 
For the solution of this model see e.g. Appendix A of ref. [3] or textbooks; 
e.g. ref. [T7] . 

In the spin-wave approximation, vortices that drive the KT phase tran- 
sition are absent. A careful analysis shows that they are, in an RG-sense, 
irrelevant for I3sw > f OH]. I.e. the large distance behaviour of systems at 
the KT-transition is given by the spin- wave model at [3sw = f • 

For the discussion of the RG-fiow in the neighbourhood of the KT- 
transition it is convenient to define x = ir/Ssw ~ 2. At the KT-transition, x 
behaves as [3l H] 

where C is an integration constant and / is a length scale. In the case of 
finite size scaling, we identify the lattice size L with this scale. 

Leading corrections to the asymptotic behaviour of correlation functions 
at the KT-transition can be obtained by computing the correlation function 
in the spin-wave approximation for Psw given by eq. 

Along these lines we obtained [18] 



T nfi-^fi^nsiTsiQ , 0.318899454... , 

JLL2,transition = 0.63650817819... H ■ ^ „ h ... (23) 

m L + C 

for the helicity modulus. Note that in the literature (see e.g. eq. (4) of ref. 
[T9] ) often T = I + ^ (i^l^c) • • • is given. The small difference is due to the 
fact that in ref. [18] contributions from configurations with non-zero winding 
number are taken into account. For the second moment correlation length 
over the lattice size we obtained [T8] 

^2nd 



= 0.7506912... + »™- +... . (24) 

L'^, transition lu L + G 
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The subscript refers to the fact that the result holds for a lattice of the 
size Li = L2 = L with periodic boundary conditions at the KT-transition. 
Furthermore the couplings between the spins have to be the same in both 
directions as it is the case here. 

Computing four-point functions {{sx ■ Sy){su ■ for all y — x, u — x and 
V — X from the propagator of the Gaussian model requires a triple sum over 
all points of the lattice. Therefore we could not reach sufficiently large lattice 
sizes this way. To avoid this problem, we performed Monte Carlo simulations 
of the Gaussian model instead. 

4.1 Monte Carlo Simulation of the Gaussian Model 

First note that Psw in eq. fl2T|) can be absorbed into the field variable by 

1/2 

(f)x = (^gyyOx- In terms of the new field variable the Hamiltonian becomes 



where ipxi^i, ^2) = 27r(?T,iXi +n2X2) / L takes into account windings of the spin 
configuration, where rii and ^2 are integers. As discussed in section 3.2 of ref. 
|18j . expectation values of an observable A in the spin wave approximation 
on a lattice with periodic boundary conditions are given by 




(25) 



The spin configurations are given by 



Sx{ni,n2) = (cos{(3g^'^(f)x + 'ipx{ni,n2)),sm{(3s^^(f)x + tlJx{ni,n2))) , (26) 



En„n.^K>^2)(A({^ln„nJ) 



(27) 



where 



{A{{s}n„n,)) 



I Uxl^^MEx 0.) expi-Hswm)) 



(28) 



The weights for the different winding numbers are given by 



W{ni, ^2) = exp [-2TT^Psw[nl + nj 



(29) 
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Note that these weights do not depend on L and therefore also non-zero 
winding numbers contribute to the asymptotic behaviour. In our simu- 
lations we have only taken into account the winding numbers {ni,n2) = 
(— 1, 0), (1, 0), (0, — 1) and (0,1). The weight for these winding numbers is 
0.000003487. . . at Psw = S/vr. This weight is rather small. However it turns 
out that at the level of accuracy that we have reached these windings have 
to be taken into account. Higher winding numbers can be safely ignored. 

In momentum space, the degrees of freedom of the Gaussian model de- 
couple. Therefore one can directly generate a field configuration, using e.g. 
the Box-Muller algorithm, in momentum space with the correct probabil- 
ity density. Then one can perform a Fourier transformation to obtain the 
configuration in real space. 

Since we had the program code available, we used a different approach: 
We generated the configurations directly in real space, using a mixture of the 
local Metropolis, the overrelaxation and a single cluster version [20] of the 
valley to mountain reflection (VMR) algorithm [21]. Due to the use of the 
cluster algorithm, there should be no slowing down and therefore our choice 
should have a similar performance as the one sketched above. 

In an elementary step of the Metropolis algorithm we propose to change 
the field at the site x as 

0L = 0, + 2(r-l/2), (30) 

were r is a random number that is uniformly distributed in [0, 1]. The pro- 
posal is accepted with the probability A = min[l, exp(— Aif)]. An elemen- 
tary step of the overrelaxation algorithm is given by 

<Px = 2 ^ (j)x , (31) 

where y.nn.x indicates that y is a nearest neighbour of x. Note that this 
update does not change the value of the Hamiltonian. In both cases, the 
Metropolis update and the overrelaxation update, we go through the lattice 
in lexicographic order. Going through the lattice once with the local update 
is called sweep. 

The variant of the VMR algorithm that we have used here is given by the 
following steps: 

• Chose randomly a site x of the lattice. The reference hight for the 
update is then given by Hq = (p^- All fields 0^ that reside on sites z 
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that belong to the cluster are updated as 

0; = 2ho - 0. . (32) 

I.e. the field is refiected at Hq. 

• Chose randomly a site y of the lattice as the starting point of the single 
cluster. The cluster consists of all sites z that are connected with y by 
a chain of frozen links. The links that are not frozen are called deleted. 
The probability to delete a link f , /i is given by 

Pd{v,fi) = min[l,exp(-2(0„ - /io)(0^+a - ^0))] • (33) 

Note that in this equation both and (pv+ji are taken before the up- 
date. 



We performed the updates in the following sequence: One Metropolis 
sweep, one overrelaxation sweep and finally a VMR single cluster update. 
The average size of the single cluster is about 1/3 of the lattice, independent 
of the lattice size. Integrated autocorrelation times of the magnetic suscepti- 
bility ( fT3f2711 at Psw = 2/11 are about 2.8 for L = 32 and increase to about 5 
for L = 512, where the time unit is one update sequence as specified above. 

In our production runs, the measurements are separated by 5 such se- 
quences. We performed 20, 20, 12, 5 and 2.2 xlO'' measurements for L = 32, 
64, 128, 256 and 512, respectively. In total, these simulations took about 10 
month of CPU time on a 3 GHz Pentium 4 CPU. 

In order to compute lim^^^oo at the KT-transition, we have set Psw = 

^ in eq. (|26ll . In addition we have used f3sw = + O.OOOlj and (3sw = 

(^^J^ ~ 0.0001 j to compute the derivative of limi^oo U4 with respect to (3sw 

at jSsw = f • Note that we have used the same set of configurations for the 
three values of Psw- Our results are summarized in tabled! 

A brief remark on non-zero winding numbers: Taking only zero winding 
configurations, f/4 is about 0.000013 smaller than with the non-zero wind- 
ing numbers taken into account. I.e. the contribution of non-zero winding 
numbers is a little larger than our statistical error. 

In order to obtain a result for the limit L ^ 00 we have fitted the results 
to the ansatz X{L) = X(oo) + cL~^, where X is either the Binder cumulant 
or its derivative. Our final results are taken from the fit that includes the 
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Table 1 : Results for the Binder cumulant U4 in the spin wave approximation 
at Psw = 2/7r and its derivative with respect to Psw- the last row we give 
the result for the limit L ^ 00. 



L 


u. 


dU,/dPsw 


32 


1.018554(2) 


-0.057529(7) 


64 


1.018298(2) 


-0.056679(7) 


128 


1.018217(3) 


-0.056394(9) 


256 


1.018200(5) 


-0.056332(15) 


512 


1.018199(8) 


-0.056297(25) 


00 


1.018192(6) 


-0.056303(16) 



lattice sizes L = 128, 256 and 512. These results are consistent within error 
bars with those obtained from L = 64 and 128. Hence the systematic error 
due to higher order corrections should be quite small. 

Plugging the result for the derivative of f/4 into eq. (l22l) we arrive at 

U,,LMransiUon = 1.018192(6) - "'"^^^^^^ + - , (34) 

where we should note again that the result only holds for a lattice with 
Li = L2 = L with periodic boundary conditions at the KT-transition. Fur- 
thermore the couplings between the spins have to be the same in both direc- 
tions as it is the case here. 



5 Monte Carlo Simulations 

In this section we discuss the details of our Monte Carlo simulations of the 
XY, the Villain and the dual of the ASOS model. 

5.1 Details of the Simulations 

We have simulated the XY model, the Villain and the dual of the ASOS 
model at the best estimates of the KT-temperature given in eq. (fTTl) . In 
our simulations we have used the single cluster algorithm [20]. Let us briefly 
summarize the steps of the cluster update: 
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• Chose a direction: 



= (cos(27rr), sin(27rr)) , (35) 

where r is a random number which is uniformly distributed in [0, 1] 0. 
All spins Sx that reside on sites x that belong to the cluster will be 
updated as 

Sx = Sx - 2{d ■ Sx)d , (36) 

Chose randomly a site y as starting point of the single cluster. All 
sites X belong to the cluster that are connected by a frozen chain of 
links with the site y. A link that is not frozen is called deleted. The 
probability to delete a link x, /x is given by 



Pd{x,^) = min 



w [sx- Sx+f, - 2{d ■ Sx){d ■ Sx+^) 

1,^ . 

wis., 



(37) 



Note that in this equation both Sx and Sx+^ are taken before the update. 

Our first set of simulations of the XY-model was performed with the single 
cluster algorithm alone. In these simulations we have used the G05CAF 
random number generator of the NAG-library. We simulated the lattice 
sizes L = 16, 32, 64, 128, 256, 512, 1024, 2048 and 4096. For L < 1024 we 
performed 5 x 10^ measurements. For L = 2048 and L = 4096 we performed 
2.5 X 10^ and 1 x 10^ measurements, respectively. In all these cases, we 
performed 10 single cluster updates for one measurement. In total, these 
simulations took a bit more than 7 month of CPU time on a 3 GHz Pentium 
4 CPU. 

Later we added simulations for L = 12, 24, 48, 96, 192, 384 and 768. 
As random number generator we have used here the SIMD-oriented Fast 
Mersenne Twister (SFMT) [22] generator. In particular we use the function 
genrand^es53() that produces double precision output. 

In these simulations we performed overrelaxation updates in addition to 
the cluster updates. An overrelaxation update of the spin at the site x is 
given by 

sj = Sx- 2^^^^Sx , (38) 



■^In analogy with the VMR algorithm for the Gaussian model, one could chose randomly 
some site z of the lattice and then chose d — s^. We did note test this alternative. 
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where 

y.nn.x 

where y.nn.x means that y is a nearest neighbour of x. It is easy to check 
that this update keeps the value of the Hamiltonian constant. 

For L < 192 we performed 10^ measurements. For each measurement 
we performed 10 single cluster updates and 10 overrelaxation sweeps. For 
L = 384 and 768 we performed 9 x 10^ measurements. In these two cases 
we performed 10 single cluster updates and 5 overrelaxation sweeps for each 
measurement. In total, these simulations took a bit more than 4 month of 
CPU time on a 3 GHz Pentium 4 CPU. 

In the case of the Villain model and the dual of the ASOS model, we 
have implemented the weight function w{sx ■ Sx+fi) for the links as a table 
with 100001 entries for the arguments -1, -0.99998, -0.99996, 1. Then, 
for a given value of the argument, we linearly interpolate between the two 
closest values contained in the table. This way, the maximal error in ratios 
of the weight for links is about 10~^ and 1.4 x 10~^ for the Villain model and 
the dual of the ASOS model, respectively. For all simulations of the Villain 
model and the dual of the ASOS model we have used the SFMT random 
number generator. 

The Villain model was simulated on lattices of the size L — 16,32,64,128, 
256,512,1024 and 2048. For L < 128 we performed 10^ measurement for 
each lattice size. The statistics for the larger lattice sizes is about 8.2 x 10^, 
5.1 X 10^ 3.1 X 10^ and 1.4 x 10^ measurements for L = 256, 512, 1024 and 
2048, respectively. In all cases, we performed 10 single cluster updates per 
measurement. In total, the simulations of the Villain model took about 3 
month of CPU time on a 3 GHz Pentium 4 CPU. 

In the case of the dual of the ASOS model we have simulated lattices of 
the size L = 16,24,32,48,64,96,128,192,256,384,512,768 and 1024. Also here 
we performed 10 single cluster updates per measurement. The number of 
measurements was 2 x 10^ for all lattice sizes up to L = 512. We performed 
1.27x 10''' and 1.1 x 10''' measurements for L = 768 and L = 1024, respectively. 
In total, the simulations of the dual of the ASOS model took a little less than 
10 month of CPU time on a 3 GHz Pentium 4 CPU. 

For all cases, the integrated autocorrelation times for e.g. the magnetic 
susceptibility are close to one in units of measurements. We have discarded 
at least the first 1000 measurements for equilibration. Given the small auto- 
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Table 2: Fits of the helicity modulus T for the XY model with the ansatz lfW\) . 
Lmin Is the minimal lattice size that is included in the fit. C and A/3 are the 
parameters of the fit. For a discussion see the text. 



f-'tnin 


C 


A/5 


xVd.o.f. 


96 


1.001(6) 


0.00049(2) 


9.83 


128 


1.044(10) 


0.00038(3) 


6.22 


192 


1.067(11) 


0.00033(3) 


4.51 


256 


1.124(19) 


0.00021(4) 


2.53 


384 


1.167(24) 


0.00013(5) 


0.56 


512 


1.170(37) 


0.00013(7) 


0.74 


784 


1.191(51) 


0.00009(9) 


0.93 



correlation times this is by far sufficient. 
5.2 Analysis of the Data 

First we have analysed the data for the helicity modulus of the XY model. 
We compared these data with the behaviour given by eq. (123|) . To this 
end, we have approximated the helicity modulus in the neighbourhood of 
the simulation point by the Taylor expansion around the simulation point to 
linear order. I.e. we have used the ansatz 



T(1.1199)+1I 



0.318899454 , , 

A/5 = 0.63650817819 + — — — (40) 

/3=i.ii99 lnL + 6 



with the two parameters A/5 and C to fit our data for the helicity modulus. 
The results of these fits are summarized in table [2J In these fits we have 
used all data starting from some minimal lattice size up to the max- 
imal size available. In order to get a x^/d.o.f. close to one, Lmin = 384 is 
needed. The result for A/5 is decreasing with increasing Lmin- Starting from 
Lmin = 384, the result for A/5 is consistent within error-bars with A/5 = 0; 
i.e. our previous estimate for f^xr- Since corrections to the ansatz (HOl) are 
expected to decay slowly with increasing lattice size, it is difficult to give a 
reliable estimate of systematic errors. Therefore we abstain from giving a 
final estimate for Pkt obtained from these fits. 
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Table 3: Fits of the correlation length over the lattice size ^/L for the XY 
model using the ansatz ( Blj) . L^m is the minimal lattice size that is included 
in the fit. C and A(3 are the parameters of the fit. For a discussion see the 
text. 



Lmin 


C 


AP 


xVd.o.f. 


32 


1.63(2) 


0.00007(5) 


1.61 


48 


1.64(2) 


0.00005(6) 


1.71 


64 


1.64(3) 


0.00004(7) 


1.89 


96 


1.62(3) 


0.00009(7) 


1.82 


128 


1.64(4) 


0.00005(9) 


2.02 


192 


1.60(5) 


0.00012(10) 


1.94 



Next we have studied the second moment correlation length ^/L of the 
XY-model using the ansatz 



e(i.ii99) 1 de 

L L dp 



= 0.7506912+^. (41) 

/3=i.ii99 InL + G 



Our results are summarized in table [31 In contrast to the helicity modulus, 
the x^/d.o.f. is rather small already for Lmin = 32. Also the dependence 
of the result for A/? on L^in is much smaller than for the helicity modulus. 
Hence the amplitude of corrections to the ansatz should be much smaller than 
for the helicity modulus. On the other hand, for a given Lmin the statistical 
error of A/3 obtained from C,/L is about three times larger than that from T. 

Since A/5 as a function of Lmin behaves quite differently for ^/L and T 
we might take the difference as an estimate of the error due to the higher 
order contributions missing in the ansatze (14011411) . Starting from Lmin = 384 
the result for A/3 from helicity modulus is consistent within statistical errors 
with that obtained from ^/L for Lmin > 32. Hence the systematic error 
should not be larger than the statistical error for these Lmin- Therefore we 
quote Pkt = 1.1200(1) as final result for the inverse of the KT-transition 
temperature of the XY model from the combined analysis of T and C,/L. 

Next we have fitted ^/L for the dual of the ASOS model with the analogue 
of the ansatz fHTj) : I.e. 1.1199 is replaced by 0.80608 as argument of ^ and 
its derivative. The results are summarized in table HI The x^/d.o.f. becomes 
smaller than one starting from Lmin = 48. Also starting from L^in = 48, 
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Table 4: Analogue of table\^for the ASOS model. 





C 


AP 


xVd.o.f. 


32 


0.493(6) 


-0.000098(11) 


10.50 


48 


0.545(8) 


-0.000035(13) 


0.59 


64 


0.548(11) 


-0.000031(15) 


0.64 


96 


0.548(16) 


-0.000032(19) 


0.75 


128 


0.543(20) 


-0.000036(22) 


0.88 



Table 5: Analogue of table\^for the Villain model. 





C 


A/? 


xVd.o.f. 


16 


2.78(2) 


-0.00009(8) 


2.44 


32 


2.72(3) 


0.00006(9) 


1.30 


64 


2.69(5) 


0.00014(12) 


1.38 


128 


2.69(7) 


0.00014(17) 


1.85 



the result for A/5 depends very little on L^m- These fits suggest the new 
estimate (3kt = 0.80605(2). This result is consistent with that of ref . [Ilj . 

Finally we have fitted the data for ^/L for the Villain model with the 
analogue of the ansatz (I^Tj) . The results are summarized in table |3 We 
get x^/d.o.f. close to one starting from Lmin = 32. These fits suggest as 
estimate for the inverse of the KT-transition temperature Pkt = 0.7517(2). 
Again this result is consistent with that of ref. [13] . 

Next let us discuss the results for the parameter C in eq. (HTi) . As our final 
results we take Cxy = 1.62(3) from the fit with Lmm = 96 for the XY model, 
Cviiiain = 2.69(5) from Lmin = 64 for the Villain model and Casos = 0.54(2) 
from Lmin = 128 for the dual of the ASOS model. The differences of these 
constants for different models should be given by the logarithm of ratios of 
the scale factors bm summarized in eq. ffTTl) . In particular we get Cvaiain — 
Casos = 2.15(5) to be compared with log{bm,Asos/bm,viiiain) = 2.16(3), 
Cviiiain-CxY = 1.07(6) to be compared with log(bm,xY/bmyiiiain) = 1.07(3) 
and finally Cxy-Casos = 1-08(4) to be compared with \og(bm,Asos/bm,XY) = 
1.10(2). I.e. the results for C are fully consistent with the relative scale fac- 
tors bm obtained in ref.|14j. 
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Next we discuss the behaviour of the Binder cumulant. In fig. [T] we 
have plotted the Binder cumulant at the old estimates of (3kt summarized 
in eq. ([H]) as a function of Lxy = Lmodeibm,modei/bm,XY- The data from the 
three different models fall nicely on top of each other, indicating a universal 
behaviour. From the leading behaviour, eq. flMl) . we would expect that the 
Binder cumulant is increasing with increasing lattice size L. However we ob- 
serve quite the opposite. For all three models, it is decreasing and apparently 
reaches a stable value ~ 1.0170 for large lattice sizes. 

Within the KT-picture this might be explained by large 0(l/(lnL)^) 
correction due to the effect of vortex pairs with a distance ^ L/2. Configu- 
rations with such a vortex pair have a much smaller magnetisation than those 
without. Therefore the appearance of such vortices enlarges the value of the 
Binder cumulant U4. Since /Ssw is monotonically decreasing with increasing 
L, the effect seen in the Binder cumulant can only be explained by the effect 
of vortices. 

Based on this observation, we have fitted the Binder cumulant at the 
estimate of Pkt given in eq. ffTT]) with the ansatz 

where C and C2 are the free parameters of the fit. Results for these fits for 
the XY model, the Villain model and the dual of the ASOS model are given 
in tables [6l [7] and [HI respectively. 

X^/d.o.f. becomes smaller close to one for Lmin = 96, 32 and 96 for the 
XY model, the Villain model and the dual of the ASOS model, respectively. 
For these values of Lmin, the results for C2 of the XY model and the Villain 
model are fully consistent, supporting the universality of the value of this 
coefficient. The result for the dual of the ASOS model is slightly smaller 
than that for the other two models. This small difference might be caused by 
higher order corrections or might be due to the uncertainty in the estimate of 
(3kt- To check the latter, we have repeated the fits for the central estimate 
for PxT given in eq. ffTT]) plus or minus the error-bar. Let us give only a 
few results with x^/d.o.f. smaller than two: For the XY at (3 = 1.1200 and 
Lmin = 96 we get C = 0.006(36) and ca = 0.06764(31). For the Villain model 
at p = 0.75177 and L^in = 32 we get C = 1.04o]35) and C2 = 0.06699(32). 
And finally, for the dual of the ASOS model at j3 = 0.80606 and L^in = 96 
we get C = -1.193(10) and Ca = 0.06730(17). I.e. the small difference in 
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Table 6: Fits of the Binder cumulant of the XY model at (5 = 1.1199 with 
the ansatz 



Lmin 


C 


C2 


xV d.o.f. 


32 


-0.002(14) 


0.06821(21) 


2.55 


48 


-0.024(18) 


0.06801(24) 


2.47 


64 


-0.012(27) 


0.06810(28) 


2.71 


96 


0.069(37) 


0.06858(32) 


1.78 


128 


0.095(39) 


0.06875(39) 


1.95 



Table 7: Fits of the Binder cumulant of the Villain model at (3 = 0.75154 
with the ansatz 



Lmin 


C 


C2 


XV d.o.f. 


16 


1.030(19) 


0.06794(26) 


3.32 


32 


1.146(37) 


0.06874(34) 


1.12 


64 


1.174(72) 


0.06886(43) 


1.35 



the values of C2 observed above can indeed be explained by the uncertainty 
of the estimate of Pkt- 

Finally let us compare the differences in C for the three models with the 
results for bm given in eq. ffTTl) . Using the results obtained for Lmin = 96, 32 
and 96 for the XY-model at p = 1.1199, the Villain model at p = 0.75154 and 
the dual of the ASOS model at P = 0.80608, respectively, we get: Cvuiain — 
Casos = 2.33(4), CxY - Casos = 1-25(4) and Cvuiain - Cxy = 1-08(5). 
The two differences that involve the dual of the ASOS model are by about 
0.15 too large compared with the result obtained from bm- The difference 
Cvuiain ~ CxY IS fuUy cousisteut with that obtained from bm- Taking the 
result for the XY-model at p = 1.1200, the Villain model at p = 0.75177 and 
the dual of the ASOS model at P = 0.80606 instead, we get Cvuiain— Casos = 
2.23(4), Cxy - Casos = 1-20(4) and Cv^lla^n - Cxy = 1-03(5). I.e. the 
deviation from the expected result is much reduced. 

Finally, let us discuss the behaviour of the Binder cumulant at the KT- 
transition. Taking e.g. the result of the fit for the dual of the ASOS model 
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10000 
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1.021 
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1.018 
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10000 



Figure 1: Data for the dual of the ASOS model, the XY model and the 
Villain model at the KT transition. (See eqs. ( ITTl) ) In the upper figure we 
plot S,/L as a function of the lattice size Lxy of the XY model. Following 
eqs. f fTTl) . the lattice size of the dual of the ASOS model and the Villain 
model are rescaled as Lxy = 0.3345 x Lasos and Lxy = 2.906 x Lmiiain- 
The solid lines give ^/L and U4, obtained from eqs. fH6l) . 
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Table 8: Fits of the Binder cumulant of the dual of the ASOS model at 
p = 0.80608 with the ansatz (ig). 



Lmin 


C 


C2 


xV d.o.f. 


32 


-0.933(3) 


0.07197(10) 


164.8 


48 


-1.069(5) 


0.06936(12) 


24.9 


64 


-1.127(7) 


0.06844(13) 


7.17 


96 


-1.180(10) 


0.06769(17) 


0.67 


128 


-1.184(14) 


0.06764(21) 


0.77 



for Lmin = 96 and (3 = 0.80608 one obtains that the minimum of the Binder 
cumulant at the transition is U^^min ~ 1.017006 at L ^ 6200. To see the 
Binder cumulant increase again, rather large lattice sizes are needed: To get 
a value 1.01701 we have to go up to L ^ 11000 and 1.01705 is reached at 
about L = 38000. I.e. it is impossible with todays computing resources to 
see explicitly from Monte Carlo simulations that the Binder cumulant at the 
KT-transition is an increasing function of L for sufficiently large L. 

6 The matching method 

The RG-flow in the neighbourhood of the KT-transition is given, up to ir- 
relevant scaling fields, in terms of two coupling constants, x = nPsw ~ 2 
and the fugacity 2; [3, H]. Therefore, in finite size scaling we can map the 
KT-fiow using two different phenomenological couplings Ri{x{L), z{L)) and 
R2{x{L), z{L)). In the ideal case, 

Ri{x,z) = Ri{0,0) +cix + 0{z^,xz) (43) 
R2ix,z) = R2i0,0)+C2Z + O{x^,xz) . (44) 

In the SOS representation, such pairs of phenomenological couplings can 
be easily found [121 [Hj. In the XY representation, on a finite lattice with 
periodic boundary conditions vortices come in pairs. Therefore only phe- 
nomenological couplings of the type (H3|) can be constructed. H 

§We have not studied free boundary conditions, where also a single vortex might exist. 
For such boundary conditions we expect larger power law corrections than for periodic 
boundary conditions. 
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Since both the hehcity modulus T and the second moment correlation 
length over the lattice size ^/L can be well fitted by eqs. fl23f24l) it is very 
hard numerically to observe any explicit ^-dependence of these quantities. 
Therefore it is actually fortunate that the Binder cumulant displays a 
contribution with an apparently large amplitude. 

Therefore we can at least map the RG-flow in an intermediate regime 
of lattice sizes using ^/L or T as first phenomenological coupling and as 
second phenomenological coupling. Here, an intermediate regime of lattice 
sizes means that L should be large enough such that power law corrections 
due to irrelevant scaling fields can be ignored and on the other hand L is 
still small enough such that ^^-contributions in the numerical data of are 
much larger than the statistical error. 

The purpose of this mapping is to allow for an accurate determination 
of (3kt from relatively small lattices. This can be achieved by solving the 
following set of equations: 



where NM stands for new model, where (3kt,nm is not known and SM stands 
for solved model. In the case of the solved model, I3kt should be at least 
determined numerically to a high precision and also the phenomenological 
couplings should be determined accurately. The unknowns of the eqs. f H5i) 
are Pkt,nm of the new model and the matching factor 6^- 

The right sides of eqs. fHSl) might be given by the results of our simulations 
discussed in the previous section. These results could be represented by 
tables. To obtain values of the phenomenological couplings for any lattice 
size Ljnin < L < Lmax, where Lmin and L^ax are the minimal and maximal 
lattice size that are simulated, one might e.g. interpolate the values contained 
in the table in linearly ln(L). 

Here instead we follow an approach that is simpler to implement. As we 
have seen in the previous section, our numerical data for ^/L and 1/4^ are 
well described by the ansatze (14 If 42 1) . Therefore we suggest to use these, 
along with the results for the parameters, as right side of eqs. Since in 

ref. |T3j the transition temperature is most accurately determined for the 
dual of the ASOS model, we use the data of this model for the matching. 



Ri,nm{L, Pkt,nm) 

R2 



Rl,SM{bmL, PkT,Sm) 
R2,SM{bmL, PkT,Sm) 



(45) 
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In particular we take the result obtained from L^m = 96 at 
(3 = 0.80608 : 

.nio.no 0.017922 0.06769 

U4,AsosiL) = 1.018192 -— — - + 



lnL-1.18 (InL- 1.18)2 

^2ndAsos{L) 0.212430 , , 

^2nd,AsubK J ^ 0.7506912 +— — . 46) 

L lnL + 0.573 ^ ^ 

For consistency we also take C,/L at P = 0.80608 and Lmin = 96. Note that 
in the present context, it is only important that the equation describes the 
data well for a certain range of lattice sizes (96 < L < 1024), where the 
matching takes place. 

To check the error due to the uncertainty of Pkt one should repeat the 
matching with the results obtained at 

P = 0.80606 : 

. n.o.no 0.017922 0.06730 
U4,AsosiL) = 1.018192 -— + 



InL- 1.193 (InL- 1.193)2 

^^-^'^^osiL) ^ 0.7506912 + ^:^1^. (47) 
L In L + 0.557 ^ ^ 

Here we do not quote statistical errors and covariances of C for ^/L and 
C2 and C for U4. Note that applying the matching with eqs. fl46f47p it makes 
no sense to produce larger statistics for the new model than the one for the 
ASOS model here. Staying well below our present statistics, say 10^ or less 
measurements for the new model, the statistical error of the coefficients in 
eqs. ( H^flTl) can be safely ignored. Since we are aiming at models that are 
harder to simulate than the models discussed in the present work, this is no 
serious limitation. 



7 A first application: the (p^ model at A = 2.1 



We have tested the new method at the example of the 2-component 0^ model 
on the square lattice. The classical Hamiltonian is given by 



(48) 



x,n 
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where the field variable (p^ is a vector with two real components. In our 
convention, the Boltzmann factor is given by exp(— if^4). Since l^^^l is not 
restricted to one, the model can not be mapped into an SOS model. For our 
test we have chosen the value A = 2.1, which is a good approximation of the 
improved value of the three dimensional model [23] . In the two dimensional 
case this value has no particular meaning. 

We have simulated the model with a mixture of the single cluster algo- 
rithm, the overrelaxation update and a local Metropolis update. Here the 
local Metropolis update is needed to change \(t)x\- The single cluster algo- 
rithm, and the overrelaxation update have exactly the same form as in the 
case of the XY model. The proposal for the Metropolis update is given by 

4>f' = (49) 

with s = 3. r is a random number with a uniform distribution in the interval 
[0, 1]. is the sum over the nearest neighbours of x. The proposal for the 
second component of the field is constructed such that H^4 is only changed by 
little. The proposal is accepted with the probability A = min[l, exp{—AH)]. 
This Metropolis update is followed by a second update at the same site, where 
the role of the two components of the field is exchanged. These updates 
were performed the following sequence: One Metropolis sweep, one sweep 
with the overrelaxation algorithm, 10 single cluster updates, followed by 5 
sweeps with the overrelaxation algorithm. After such a sequence of updates, 
a measurement of the observables is performed. 

We have simulated L = 32 at (3 = 1.075 and 1.08, L = 6A at (3 = 1.078 
and 1.08, L = 128 at /5 = 1.077 and P = 1.079 and L = 256 at /3 = 1.0775 and 
f3 = 1.0785. In all cases we have performed 10^ measurements. For fitting 
the data with the ansatze (14011411) and for the matching (H5|) . we linearly 
interpolated the values of T, ^/L and f/4 in /?. 

First we have fitted the helicity modulus to the ansatz (1401) . where (3kt 
and the integration constant C are the free parameters. The results of these 
fits are summarized in table [91 x^/d-o.f. is quite large and the result for Pkt 
is decreasing with increasing Lmin- 

Next we have fitted the data for C,/L to the ansatz fHTl) . The results 
are summarized in table [TOl The x^/d.o.f. is already smaller than two for 
Lmin = 32. But on the other hand also the statistical error of i3kt obtained 



22 



Table 9: Fitting T of the 0^ model at X = 2.1 with the analogue of the 
ansatz ( BOj) . 



Lmin 




C 


xVd.o.f. 


32 


1.07865(1) 


0.278(2) 


75.9 


64 


1.07838(2) 


0.363(7) 


22.6 


128 


1.07818(2) 


0.447(8) 





Table 10: Fitting ^/L of the (j)^ model at X = 2.1 with the analogue of the 
ansatz ( Blj) . 







C 


xVd.o.f. 


32 


1.07721(6) 


1.50(2) 


1.67 


64 


1.07739(12) 


1.41(6) 


1.39 


128 


1.07761(19) 


1.30(9) 





from ^/L is about 6 times larger than for that from the helicity modulus. In 
the case ^/L the estimate of Pkt is an increasing function of Lmm- Therefore 
we might regard the result from ^/L and Lmm = 128 as lower bound and 
that from Lmin = 128 from the helicity modulus as upper bound. Hence 
we quote Pkt = 1.0779(5) as our final estimate obtained from the combined 
analysis of T and C,/L. 

Next we applied the matching method discussed in the previous section. 

Table 11: Matching the (j)^ model at X = 2.1 with the dual of the ASOS 
model at 13 = 0.80608, eqs. g^). 



L 


Lasos 




32 


55.81(12) 


1.07892(4) 


64 


113.6(5) 


1.07822(3) 


128 


234.2(1.4) 


1.07786(2) 


256 


479.(5.) 


1.07776(2) 
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Table 12: Matching the 0*^ model at X = 2.1 with the dual of the ASOS 
model at (3 = 0.80606, eqs. g^). 



L 


Lasos 


Pkt 


32 


56.03(12) 


1.07898(4) 


64 


113.8(5) 


1.07827(3) 


128 


234.0(1.4) 


1.07790(2) 


256 


478.(5.) 


1.07780(2) 



This way we obtain an estimate of (5kt and the corresponding lattice size 
of the ASOS model from a single lattice size L. Our results for matching 
with [3 ASOS = 0.80608, eqs.fH6l). are given in table fTTl and those for Pasos = 
0.80606, eqs.(H71). in table [121 The estimates for Pkt apparently converge as 
L increases. Taking into account the difference between L = 128 and L = 256 
as well as the difference between Pasos = 0.80608 and Pasos = 0.80606 we 
arrive at our final estimate PkT = 1.0778(2). This estimate is clearly more 
precise than that obtained from fitting T and ^/L with the ansatze fl40|4ip . 

8 Summary and Conclusions 

We have studied the XY model, the Villain model and the dual of the abso- 
lute value solid-on-solid model on a square lattice at the Kosterlitz-Thouless 
transition. We have focused on the behaviour of phenomenological couplings 
like the helicity modulus T, the second moment correlation length over the 
lattice size ^/L and in particular the Binder cumulant f/4. 

Using Monte Carlo simulations of the Gaussian model on the square lat- 
tice, we have determined the asymptotic value of f/4 and the leading loga- 
rithmic correction at the Kosterlitz-Thouless transition. See eq. (1341) . 

We confronted the predictions fl23ll24ll34p for the phenomenological cou- 
plings derived from the Kosterlitz-Thouless theory with the results of our 
Monte Carlo simulations. We have reached lattice sizes of L = 4096, 2048 
and 1024 for the XY model, the Villain model and the dual of the ASOS 
model, respectively. We have generated several 10^ statistically independent 
configurations for each lattice size. 

We find that the data for the helicity modulus T and in particular the 
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second moment correlation length over the lattice size ^/L follow quite well 
the predictions fESll^Ij) . 

In contrast, the value of the Binder cumulant at the KT-transition is 
increasing with increasing lattice size, while eq. flM|) predicts that the Binder 
cumulant should decrease for sufficiently large lattice sizes. This discrepancy 
is explained by the presence of vortex pairs that have a distance of order L/2. 
This effect should be proportional to the fugacity squared. Hence it should 
lead to a subleading correction. Indeed adding a term C2/(lnL + CY to the 
ansatz allows to fit our data for the Binder cumulant for all three models 
studied. The values for C2 obtained for the three models are consistent, 
confirming the universal character of the subleading correction. From the 
results of these fits one can easily see that with the present computer resources 
it is not possible to reach sufficiently large lattices with high statistics to see 
explicitly the asymptotic behaviour of the Binder cumulant (1341) . 

The fact that behaves quite differently from ^/L and T allows us to 
set up a matching method similar to that of ref. [131 [E] for solid-on-solid 
models. For a detailed discussion see section [HI 

As a first test, we have applied the matching method to the two compo- 
nent 0^ model on the square lattice at A = 2.1. We find that using moder- 
ately large lattices (up to L = 256), we can determine the temperature of 
the KT-transition quite accurately. We show that the new matching method 
is superiour to fits of T or ^/L to the predictions (I23|24l) . 

In the near future we plan to apply the method to thin films of the three 
dimensional two component 0^ model and the dynamically diluted XY model, 
where the duality transformation to solid-on-solid models is not available. 
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